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Abstract 

Numerical simulations are performed to study the finite temperature phase tran- 
sition in the SU(2) Higgs model on the lattice. In the presently investigated range 
of the Higgs boson mass, below 50 GeV, the phase transition turns out to be of first 
order and its strength is rapidly decreasing with increasing Higgs boson mass. In 
order to control the systematic errors, we also perform studies of scaling violations 
and of finite volume effects. 



1 Introduction 

The masses of elementary particles in the Standard Model are generated via the Higgs mecha- 
nism by the non-zero vacuum expectation value of the scalar Higgs field. At high temperatures, 
above the scale of the vacuum expectation value, the Higgs mechanism is not operative, the 
symmetry of the vacuum gets restored 0. In fact, in the early universe, according to the big 
bang cosmology, matter first existed in the symmetry restored phase. As a consequence of ex- 
pansion and cooling, a non-zero vacuum expectation value of the scalar field was developed in 
the phase transition between the symmetric phase at high temperatures and the Higgs phase at 
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lower temperatures. The properties of this electroweak phase transition might have a substan- 
tial influence on the later history of the Universe. For instance, since the number of baryons 
is not conserved in the minimal standard model |J, the small baryon asymmetry of the Uni- 
verse could perhaps be created in non-equilibrium processes during a strong enough first order 
electroweak phase transition |3|, . This offers the possibility that the baryon asymmetry can 
be explained within the minimal standard model. The resolution of this question is therefore a 
major challenge for elementary particle physics. 

The standard calculational method for the study of the symmetry restoring electroweak 
phase transition is resummed perturbation theory |5|, ||, ||. In the Higgs phase perturbation 
theory is expected to work well for not very high Higgs boson masses, since the couplings are 
small. In the high temperature symmetric phase, however, the situation is similar to high 
temperature QCD: irreparable infrared singularities occur which prevent a quantitative control 
of graph resummation J9j . Since the calculation of physical quantities characterizing the phase 
transition requires the knowledge of both phases, there is a priory no reason why perturbation 
theory could provide a quantitative treatment of the electroweak phase transition. Indeed, the 
results of perturbation theory show bad convergence. 

For a better understanding several non-perturbative methods have also been tried. For 
simplicity, fermions and the U(l) gauge field are often omitted. This can be expected on 
general grounds to be a reasonable first approximation. In this way one is left with the SU(2) 
Higgs model describing the interaction of a four-component Higgs scalar field with the SU(2) 
gauge field. Possible non-perturbative approaches include a block spin procedure leading to 
evolution equations for average actions flOfl , the e-expansion at 4 — e spatial dimensions [T1J and, 
of course, numerical simulations. After pioneering works [O, O, recent numerical simulations 



concentrated on the understanding of the finite temperature behaviour of the SU(2) Higgs model 
at large Higgs boson masses near and above the W-boson mass [14]] . Another non-perturbative 
approach is based on dimensional reduction, studying the three-dimensional effective Higgs 



theory, which is obtained in the high-temperature limit [15, 16, 17] . A further simplification 
leads to an effective scalar theory fT8|, which has also been studied numerically in the reduced 



model W§ 



The non-perturbative investigations of the electroweak phase transition did not yet lead to a 
convincing unique picture. Therefore, we decided to perform a large scale numerical simulation 
of the symmetry restoring phase transition in the SU(2) Higgs model. We stay in the original 
four-dimensional theory without reduction. This has the advantage of keeping the number 
of bare parameters small and not introducing any further approximations beyond the lattice 
regularization. First results have been published in a recent letter [20]. Here we give a detailed 
description of the techniques used and include additional results. As it is known from previous 
studies |14 1 , for Higgs boson masses near and above the W-boson mass the numerical simulations 
in the original four- dimensional model are technically difficult. Therefore we restrict the present 
calculations to smaller Higgs boson masses below 50 GeV. Since this region of parameters of 
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the minimal standard model is already excluded by experiments, our present scope is merely 
theoretical because we would like to check the validity of some other theoretical approximation 
schemes, e.g. resummed perturbation theory. We plan to extend this investigation to heavier 
Higgs boson masses in future papers. 



1.1 Lattice action 

The lattice action of the SU(2) Higgs model is conventionally written as 



pi v 



i. 



X 



1. 



i 



+ E { oTr (fp+ip.) + A -Tr (<p+<p x ) - 1 -«^Tr (fpZ^U^) \ ■ (1) 



Here U Xfl denotes the SU(2) gauge link variable, U v \ is the product of four C/'s around a plaquette 
and (p x is a complex 2 ® 2 matrix in isospin space describing the Higgs scalar field and satisfying 

ft = T 2f T x r 2 . (2) 

The bare parameters in the action are (3 = 4:/g 2 for the gauge coupling, A for the scalar 
quartic coupling and k for the scalar hopping parameter related to the bare mass square /Xq 
by /1q = (1 — 2A)/t~ 1 — 8. Throughout this paper we set the lattice spacing to one (a = 1), 
therefore all the masses and correlation lengths etc. will always be given in lattice units, unless 
otherwise stated. 

In order to fix the physical parameters in a numerical simulation one has to define and 
compute some suitable renormalized quantities at zero temperature. The physical Higgs mass 
Mh and W-mass My/ can be extracted from correlation functions of different quantities (see 
section ||). The renormalized gauge coupling can be determined from the static potential of an 
external SU(2) charge pair, measured by Wilson loops (see section H). 

Since we are interested in the study of the symmetry restoring phase transition as a func- 
tion of temperature, we use asymmetric lattices: the small temporal extensions L t = 2,3, . . . 
represent the discretized inverse temperature L t = l/(aT). The other three (spatial) extensions 
of the lattice have to be much larger, for reaching the thermodynamical limit. In order to fix 
the physical parameters at the phase transition we have to determine the zero temperature 
renormalized parameters at the phase transition points for the L t = 2, 3, . . . lattices. A renor- 
malized gauge coupling near the physical value g\ ~ 0.5 can be obtained near (3 = 8 ^T|. As 



stated before, we would like to have lighter Higgs boson masses than studied in []14| . Hence for 
the bare quartic coupling we have chosen values near A = 0.0001 (in this paper referred to as 
low) and near A = 0.0005 (referred to as high). In the present paper the inverse temperature 
in lattice units will be restricted to L t = 2 and L t = 3. Therefore the indices of the four sets of 
numerical simulations will be 
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• 12 for low A and L t = 2; 

• 13 for low A and L t = 3; 

• h2 for high A and L t = 2; 

• h3 for high A and L t = 3. 

In the next section the numerical simulation methods will be discussed. An important tool 
for the orientation in bare parameter space will be the invariant effective potential introduced in 
section |[ Then different groups of numerical simulation results will be discussed: the location 
of the phase transition points in section || masses and correlation lengths in section |j, the 
renormalized gauge coupling and the renormalization group trajectories in section [], the latent 
heat in section [7| and finally the interface tension in section |8|. The last section is devoted to 
the discussion of results and to a summary. 

2 Monte Carlo simulation 

In this section some aspects of the applied Monte Carlo simulation techniques are discussed. 
This can be skipped by readers not interested in technical details. 

The simulations have been performed on the Alenia Quadrics computers of DESY. The 
Quadrics Q16 is a massive parallel machine with SIMD[] architecture which consists of 128 
processors (nodes). Depending on the goals and features of the respective simulation, we use 
different strategies: 

• A lattice is assigned to each node. No time is wasted for the communications between 
the nodes. Limitations of memory allow this only for small enough lattice extensions. 

• The Q16 may be switched to consist of 16 independent 2 3 tori. 

• The whole machine is arranged as a three-dimensional torus. 

In this way the lattice is distributed over 1, 8 or 128 nodes and one obtains 128, 16 or 1 
independent data sets from one run, respectively. Of course, the lattice extensions have to be 
multiples of the corresponding tori. 

The Quadrics offers 32 bit floating point arithmetics. This is sufficient for most of our 
purposes, except for building global averages, when we use a simple variant of software based 
double precision arithmetics for summation. 

1 single instruction multiple data 
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2.1 Updating 

In a Monte Carlo simulation the autocorrelation of subsequent configurations is one of the 
main problems one has to deal with. The autocorrelation is usually worse at phase transitions. 
This phenomenon is called critical slowing down. Unfortunately, the continuum limit has to be 
performed in this region of parameter space (in our case k, j3 and A). 

Due to the small A-values we use, large fluctuations of the squared Higgs field length p\ = 
\Tv{ip+Lp x ) occur in the Higgs phase and this expectation value shows the largest autocorrelation 
of all investigated quantities. In order to reduce the autocorrelation, the updating has to 
propose a big change of p x , with a high probability to accept it. This is achieved by an 
overrelaxation algorithm for p x , developed by two of us ||22|. As usual, this overrelaxation 



algorithm is non-ergodic, thus we have to mix it with some ergodic updating steps. 

A further improvement of the autocorrelation was achieved by the replacement of the ergodic 
Metropolis step considered in with a heatbath algorithm proposed to us by Burkhard Bunk 
p3[ . This algorithm works on the four real components of the y^-field. As far as we know, 
this algorithm has not appeared in a publication, therefore we would like to describe it in some 
detail. 

The idea is to divide the action in two parts, a quadratic one and the rest proportional to 
A. In a heatbath step a new (p x is proposed according to the quadratic part of the action and 
the remaining parts are taken into account by an additional accept-reject procedure. Starting 
from the lattice action ([!]), the y^-dependent part can be written in the following form: 

S((p x ) = ((f x ,o - b xfl f + . . . + (<p x , 3 - b Xt3 ) 2 + X(pl - l) 2 + const. (3) 

Here we used the notations 

3 

<£x = <fx,0 1 + i E ¥x,m T m , 

771=1 

K 4 

Ko = 7J E Tr (vt+uUxv + U x _ ^x-u) , 
4 

Km = -77 J2 Tr Vft+vUxuTm + ^>t-v V t-v^m) ■ 

In this form of the action a term proportional to p x can be shifted from the quadratic part to 
the quartic part. We introduce a factor ( x to express this freedom: 

E (v>*J ~ + A {f* ~ ^( 2A " 1 + C*)) + const.' (4) 



S(lfx) = Cx 



For a good acceptance the minima of quadratic and quartic parts in ([|) should coincide. So we 
get for ( x the equation 

3 

2X\b x \ 2 = £(( x + 2A - 1) , with: \b x \ 2 = ]T b x>J b xJ . (5) 
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This cubic equation cannot be solved fast enough in an updating. Starting from the observation 
( x = 1 for 1 63. | = 1, we split ( x in ( x = 1 + e x and get the following approximate expression for 
the optimal ( x : 

( x = 1 - 2A + 2A • \b x \ 2 + 0(e 2 x ,e x \) . (6) 

This approximation works very well in a sufficient range of \b x \. In practice the average accep- 
tance of this algorithm turned out to be larger than 98%. 

For the SU (2)- variables U xn and a x = <p x /p x we use standard overrelaxation methods [^2, 



[24]] . Because the above described heatbath algorithm for (p x offers also an ergodic update for 
the angular part a x of the scalar field <p x , we need an ergodic update for U Xfl only. For this 
purpose we use the heatbath algorithm described in [^, f26fl . 

In all updatings, the random number generator proposed and implemented by Martin 
Luscher p7| is applied. It is based on an algorithm of Marsaglia and Zaman f2"8fl . The latter 



algorithm is known by the name RCARRY, if the parameters have been chosen appropriately. 
RCARRY offers an extrem long period > 10 m , but unfortunately it owns some short range 
correlations. As it has been shown ||27| , on long range, a chaotic nature of the algorithm comes 
to light. Skipping from time to time some hundreds of numbers in the sequence the corre- 
lation is practically eliminated. Due to the skip this random number generator is relatively 
slow. In order to benefit from the parallel architecture, the random number generator has to 
be initialized independently on the nodes of the machine. 

For updating the field configurations a combination of the above described five algorithms 
is used. We choose some basic sequence of elementary updatings for the different sets of field 
variables, which is repeated periodically many times. The whole sequence, which visits every 
variable at least once but usually many times, will be called sweep. The optimization of this 
basic sequence making up a sweep is a difficult but important task, which will be discussed in 
the next subsection. 

2.2 Autocorrelations 

Let us consider the autocorrelation of a quantity Q[U Xfl , <p x ] measured on a sequence of field 
configurations. Q n is the value of Q[U Xfl ,(p x ] measured on the n-th configuration and Q = 
-k J2n=i Qn is the average over the configurations. The autocorrelation function for this quantity 
is defined as: 

i N-t 

r Q (t) = Jim — £ {Q n+ t - Q) (Qn - 0) . (7) 

By the use of the integrated autocorrelation time 

r " 1 V ^ (t) (R) 
1 t=-oo i Ql U J 

the statistical error can be estimated from the variance: 

al^2r mt>Q ^-Q 2 ) . (9) 



lattice: 16 <g> 16 <8> 16 <g> 32 
j8= 8 k= 0.12885 X= 0.0005 




10 20 30 40 



t in sweeps 

Figure 1: An example of autocorrelation function in the Higgs phase at T = 0. 
The curves for T p 2(t) and Tcf(t) coincide. 



We investigated the integrated autocorrelation time for four characteristic quantities: p 2 , 
P(x+Lt/2)Pxi U p i and the largest calculated Wilson loop. The second quantity characterizes the 
correlation function of p 2 at the largest distance. In zero temperature simulations the largest 
Wilson loop had the size L s /2 ® L t /2, with L s and L t denoting the extensions of the lattice in 
time and space direction, respectively. At finite temperature only the 1 ® 1 Wilson loop was 
considered. For Wilson loops not every orientation was taken: two sides were always in the 
direction of the largest lattice extension. Further on, we refer to the correlation function as Cf 
and to the Wilson loop as Wl. If equation (H) is evaluated on a finite sequence of configurations, 
one has to decide where to truncate the sum over t. As mentioned before, the largest r int - values 
were found for pi, so we truncated at the first zero of T p 2(t). 

On larger lattices we usually had 16 independent configurations in the computer. They were 
evaluated separately and an estimate for the statistical error of the integrated autocorrelation 
time was obtained from the variance. 

We made some investigations how to optimize the autocorrelation by changing the number 
of calls of the various updatings in the complete sweeps but we did not try to optimize the 
ordering of the algorithms. Of course, we optimized the autocorrelation in CPU time since the 
above changes affect the time requirements of the complete sweeps. It should be mentioned that 
the optimal number of calls depends on lattice size and parameter range. The measurement 
routines were called after each sweep. 
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lattice: 3 <g> 24 <g> 24 ® 96 
/y = 8.15 k= 0.1281 X = 0.00011 




10 20 30 40 
t in sweeps 
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lattice: 3 g> 24 <g> 24 96 
p= 8.15 k= 0.1281 X = 0.00011 




Higgs phase 



i i i i I 



li i 



0.2 



100 200 300 
t in sweeps 



0.1 



400 



Figure 2: Autocorrelation functions at the phase transition line on a 3 • 24 2 • 96 
lattice. The left picture refers to the symmetric phase, the right one to the Higgs 
phase. Evaluated sweeps: 32000 sweeps in the symmetric phase and 80000 sweeps 
in the Higgs phase. In both cases the same updating scheme was applied. In the 
symmetric phase T int p 2 = 21 ± 4 sweeps and in the Higgs phase T intjf> 2 = 217 ± 66 
sweeps. 



The autocorrelation function was investigated in the Higgs phase both for T = and for 
finite values of T. In the symmetric phase only finite T was considered. A typical example of 
autocorrelation in the Higgs phase is shown by fig. |l| The autocorrelation function T p 2(t) is 
to a good approximation a single exponential. There is no significant difference between Tcf(t) 
and T p 2 (t) . The autocorrelation functions for quantities depending only on U Xfl show a fast fall 
off for t values very small compared to T int p 2. For larger values of t the exponential descent of 
Tjj (t) and Fwi(t) is the same as of T p 2. We always found the initial fall off to be larger for 
T W i(t) than for T Upl (t). 

A comparison of autocorrelations in the two phases is given in fig. |[ In the symmetric 
phase at finite T we found the autocorrelation time for the quantities U p i and Wl to be less 
than 1. The left hand side of figure ^| displays the extremely fast descent of these autocorrelation 
functions. The behaviour of T p 2(t) differs from the behaviour in the Higgs phase, because there 
is a strong curvature in the logarithmic plot. This could be a signal for a dense spectrum of 
states. The integrated autocorrelation time in the symmetric phase is usually much smaller 
than the one in the Higgs phase at the same parameters and lattice extensions. 
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Table 1: Integrated autocorrelation time T inttP 2 of p\ for 5 different updating 
schemes. The lattice size is 16 3 • 32, the parameters are /3 — 8 and A = 0.0005. 
The first item is at k = 0.12885 and all the others are at k — 0.1289. For each item 
more than 32000 sweeps were evaluated. 



heatbath 


overrelaxation 


Tint,p2 


U X f_l 




U X fX 


ot x 


px 


in sec. of CPU 


in sweeps 


1 


4 


3 


3 


1 


24.2 ±0.8 


11.4 ±0.4 


3 


3 


3 


3 


1 


33.3 ±2.3 


11.4 ±0.8 


1 


1 


3 


3 


1 


34.1 ±2.3 


19.6 ± 1.3 


1 


1 


6 


6 


1 


60.5 ±6.2 


22.7 ±2.3 


1 


1 


3 


3 


3 


62.1 ±4.7 


28.0 ±2.1 



On a 16 3 ■ 32 lattice with parameters (3 = 8, k = 0.1289 and A = 0.0005 we compared four 
different compositions of the complete sweep. These parameters give a point in the Higgs phase. 
The results for the largest autocorrelation time Ti ntiP 2 are given in table |I] in sweeps and in CPU 
seconds of Q16, assuming that the lattice is distributed on the whole machine. A comparison 
of the third and fourth rows shows that more work on the SU(2) variables has no influence on 
the autocorrelation: the autocorrelation time in sweeps is about the same for both updating 
schemes. The fact that more overrelaxation for p x does not lead to a better autocorrelation is 
plausible 

The autocorrelation measured in sweeps decreases significantly if more heatbath is called 
but, due to the time needed for the heatbath algorithms, there is no significant difference 
between the second and the third row of table [l], if the autocorrelation time is measured in 
CPU time. 

The updating scheme in the first row of table [1] is the best combination we found. It was 
therefore used in many points. Because of the large autocorrelations in the Higgs phase at finite 
temperatures, which were typically about 10 times longer in sweeps than at T = 0, we could 
not compare different updating schemes there. Another comparison of updating schemes was 
performed on a 18 3 ■ 36 lattice with parameters (3 = 8.15, k = 0.1281 and A = 0.00011. The 
conclusions were very similar. 

2.3 Multicanonical simulation 

An important problem of Monte Carlo simulations of a system with first order phase tran- 
sition is the supercritical slowing down. At the transition point the tunneling rate between 
the two phases is exponentially suppressed for any local update algorithm (e.g. overrelaxation, 
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heatbath). To overcome this problem the multicanonical algorithm was developed [p9 ]. The 
basic idea is an enhancement of the mixed states, which are suppressed due to the additional 
free energy of the interfaces. This enhancement is reached by an extra term in the action, i.e. 
S — > S + f{0). This term can be a function of any order parameter O. The easiest way is to 
use the action and a continuous function f(S) = faS + ®k with constant fa, ®k for S in the 
interval I k = (S' fe ,5' fe+1 ]. In fact, instead of the lattice action in eq. ([]]), we used as an order 
parameter the modified action 

S log = S[U,<f]-3^og(p x ) , (10) 

X 

which is natural to take if p x is used as an integration variable in the path integral. This 
choice is particularly convenient, since all the overrelaxation algorithms (p x , a x , U xft ) can be 
used without changes. The intervals I k and the parameters a k and fa are chosen in such a 
way that the multicanonical probability distribution P™ c is nearly flat. This is achieved if 
f{Siog) ~ log(P_L) between the two maxima and is constant elsewhere. Here Pl is the canonical 
probability distribution of the action Si og . The distribution Pl is obtained in a multicanonical 
simulation by reweighting P™ c with exp(faSi og + a k ). 

In practice a first choice for the multicanonical parameters is made and they are optimized 
afterwards. If necessary, the procedure is repeated until P|" c becomes flat. A first guess can 
be obtained from smaller lattices. In this way the distribution of the action Si og and the link 
variable L v , defined in eq. fl24|) , was measured on 2 • 4 2 • 64 and 2 • 4 2 • 128 lattices at the "low" 
value of the quartic coupling (A = 0.0001). 

For larger lattices two problems arise. The parameters have to be tuned very precisely and 
the autocorrelation times become even for optimally tuned values very large, of the order of 
(9(10000) sweeps. To solve these problems we combined the multicanonical method with the 
constrained simulation method [[Kj. In what follows we call this way of simulation constrained- 



multicanonical method. 

We divide the interval between the two maxima of Pl into subintervals. These are chosen 
to have an overlap with their neighbours. Starting in one phase we tune the multicanonical 
parameters such that P™ c is flat in a given subinterval and suppressed elsewhere. This means 
that f(Si og ) is approximately equal to log(Pt) in this subinterval and increases rapidly beyond 
the boundaries. By moving the subinterval one is going from one maximum to the other. At 
the end every set is reweighted. In case of large overlaps between neighbouring intervals the 
absolute normalization can be obtained with small errors. 

To ensure that this method yields the same result as the pure multicanonical one, we 
performed simulations using both methods on 2 • 4 2 • 128 lattice. The results coincide within 
statistical errors. Only the constrained-multicanonical algorithm has been used for simulations 
on 2 • 8 2 • 128 lattice. 

The technical realization of the multicanonical approach by the Metropolis algorithm is 
straightforward. As it has been emphasized above, due to our special choice of the order 
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parameter, the overrelaxation algorithms can be used without changes. The modifications for 
the heatbath algorithms are more involved (see e.g. ||31|| ). A description of our implementation 
of heatbath algorithms for U xpi and (p x is given in the appendix. Since the heatbath algorithms 
are more efficient, we always used them instead of the Metropolis algorithms. 

In our simulations the acceptance rate for the multicanonical heatbath algorithms was very 
good: for the modified gauge algorithm at least 99% and for the yj-algorithm at least 96%. The 
overlaps for the neighbouring intervals were chosen to be approximately 40%. The number of 
subintervals was 5 for the 2 ■ 4 2 ■ 128 lattice and 13 for the 2 • 8 2 ■ 128 lattice. The autocorrelation 
times were on average about 500 sweeps for the constrained-multicanonical simulations. For 
the two smaller lattices we measured about 2000 and 7000 sweeps as autocorrelation times with 
the pure multicanonical algorithm. 

The easiest way to parallelize the multicanonical algorithm on the Quadrics Q16 machine 
is to simulate several lattices independently on each node. For any other implementation there 
is a need for communication between the different nodes for each updating step, since / is 
a function of the global action. Another disadvantage of partitioning the lattice would be a 
decrease in the acceptance rate due to simultaneous change of several variables. 



3 Invariant effective potential 

In the perturbative approach to the electroweak phase transition the most important quantity 
to compute is the effective potential. Interesting physical observables like latent heat, surface 
tension or masses can be extracted from it. Of course, these latter quantities can also be 
obtained from the non-perturbative approach of numerical lattice simulations by measuring 
suitable observables. Nevertheless, a direct comparison of the effective potential itself from both 
methods would obviously be desirable. On the lattice, however, the action is gauge invariant 
and so are the observables, as demanded by Elitzur's theorem. In perturbation theory the 
effective potential is calculated in different gauges. (The most popular one is the Landau 
gauge.) In order to compare this with lattice results one ought to fix the gauge on the lattice, 
a notoriously difficult task in particular for non-abelian gauge groups. 

A way out is the study of the gauge invariant effective potential that has been initiated 
recently [^, ^ f\ In this approach one considers composite gauge invariant operators in the 



standard Legendre transformation framework. The obvious advantage of this approach is that 
the potential can be evaluated perturbatively and, since it is gauge invariant, it can be directly 
compared to lattice simulations. It offers therefore a conceptually clean and directly accessible 
tool of confronting results obtained in perturbation theory with numerical data. In this paper 
we want to report about our first experiences with the gauge invariant potential. We will 
use it mainly for the determination of the transition points. We postpone a discussion of its 



2 Another possibility could be the effective potential computed from the lowest eigenvalue of the Laplace 
operator 1 34 . 
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renormalization and the extraction of physical quantities to a future publication. 

The starting point for the gauge invariant effective potential for the length square of the 
Higgs field is the free energy F(J) in the presence of a constant external source J 

e- QF ^ = J [dU][dbp]e- s+J ^^ , (11) 

with S the action eq. ([[]) and Q the lattice volume. From this the effective potential is obtained 
by a Legendre transformation 

V{f) = F( J(p 2 )) - p 2 J , (12) 

where 

t - Tj F(J) . (13) 

The perturbative evaluation of the gauge invariant effective potential is done in the standard 
loopwise semiclassical expansion. To proceed we choose the unitary gauge and set the angle 
of the Higgs field to a x = 1. Since the free energy is gauge invariant this has no effect on the 
effective potential. The action becomes 



2 



S = S a + { Pi + X [Pi - lj - K Yl Px+vPx Tt U w \ » ( 14 ) 
1 I J 

with S g the plaquette action for the gauge field defined in (TJ). 
A stationary point is obtained for 

J>1-8k-2A (15) 



and U x ^ the unit matrix as 



o J + 8k - 1 



The last equation is easily inverted for J(p 2 ) and the effective potential to tree level is 

Vtrelf) = (1 - 8K)p 2 + A(p 2 - l) 2 . (17) 

To get the one-loop effective potential we have to consider fluctuations around the stationary 
point (|T6|). The fluctuations to one-loop consist of a gauge part and a Higgs part and at this 
level no mixing appears. One obtains 

Vi-ta* = V tree + £ [I ln(P + m 2 ) + ~ ln(P + m 2 )} , (18) 



where the masses are related to the parameters in the lattice action (|14D and p by 



2^2-2 2 4 \ -2 /i n\ 

m g = o K 9 P , m (j> = -Xp (19) 

Z, Hi 



and the momenta in the lattice integrals are k 2 = X)u[2 — 2 cos(A; /J )]. 
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The solution of the effective potential given above corresponds to the broken phase of the 
SU(2)-Higgs model. In it was emphazised that there exists another stationary point which 
belongs to the symmetric phase of the model and which is given by p = 0. In this case the 
tree level potential is trivially zero and we have to start with the one-loop formula for the free 
energy 

with wig = (1 — 8k, — 2\)/k + J. To obtain the effective potential one has to solve eq. flT3| ) for 
J(p 2 ). In |33| the solution has been given for the three dimensional Higgs model in a closed 
form. A description of the effective potential in the symmetric phase has been found with a 
quite characteristic asymmetric shape. 

In our case we have to work with lattice integrals or finite lattice sums. Then the solution can 
no longer be given in a closed form. However, one can perform the Legendre transformation and 
solve eq. (T3) numerically. The result of this procedure for the points where our simulations 
are performed confirm the general shape of the potential in the symmetric phase and are 
qualitatively in agreement with the potential extracted from the distributions of p 2 values from 
the simulations. However, we do not have a quantitative understanding of the symmetric phase 
yet. We hope to come back to this question in a future publication. 

A final remark concerns the lattice simulations where the gauge invariant effective potential 
is obtained from a distribution of the operator under consideration. The potential computed 
in this way is the so-called constraint effective potential |3~5|| . In the infinite volume limit this 
potential coincides with the one defined by means of the Legendre transformation above. In 
perturbation theory both approaches differ in the treatment of the zero modes. For the one- 
loop result (fl8|) this amounts to leaving out the k = mode in the finite lattice sums for the 
constraint effective potential. 



4 Phase transition points 

A numerical simulation of the SU(2)-Higgs model should start by first determining the phase 
transition points. Physically the transition is triggered by a temperature change. Keeping all 
other parameters fixed, this could only be achieved on the lattice by asymmetric couplings, 
which is possible but cumbersome. It will become clear later (see sections ^ and |6|) that in 
the parameter range we are interested in (3 and A are fixing, to a good approximation, the 
renormalized parameters. A change of k is reflected mainly in a change of the lattice spacing 
a. Therefore if one crosses the transition at fixed j3, A by changing k, the essential change is 
in the physical temperature T = l/(aL t ). (The physical volume is assumed to be large enough 
such that its change with a 3 is not important.) Thus we are looking for the phase transition in 
the hopping parameter at k = k c , for fixed (3, A. 
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lattice: 2(g>4®4<g>32 
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Figure 3: Thermal cycle exhibiting a hysteresis in p 2 . The solid line indicates the 
values of the absolut minima from the gauge invariant effective potential, the short 
dashed line the ones of the false minima. The long dashed line only connects the 
data points to guide the eye. 



4.1 One-loop effective potential and transition points 

We found that for searching the transition point the gauge invariant effective potential can 
be very helpful. It can serve as a tool to provide quite accurate information about k c which 
helps to select the K-values where simulations are then performed. We define the transition 
point k c as the K-value where the symmetric and broken minima of the gauge invariant effective 
potential are degenerate. For the computation we used the one-loop formula eq. (p~8|) for the 
broken phase and the trivial, p 2 = 0, minimum for the symmetric phase. Although this is 
certainly not the exact value for the symmetric minimum, we will see in the following that for 
small A the transition k's are in very good agreement with numerical data. 

For the computation of the gauge invariant effective potential on a finite lattice of size 
L x ■ L y ■ L z ■ L t the integrals in fll8|) have been replaced by the corresponding lattice sums. 
Following the experience in QCD ||36|| , we also used the mean field improved gauge coupling 
9 ~^ 9l\fUph with Upi the measured plaquette value. This choice for the gauge coupling 
appeared to be very useful to achieve better agreement with simulation results. 

At the values of the Higgs mass where our simulations are performed the electroweak phase- 
transition is certainly of first order. The effective potential is such that besides the absolute 
minimum a secondary minimum exists. In the simulations this phenomenon leads to hystere- 
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sis effects when the system gets stuck in the wrong minimum. Therefore hysteresis effects in 
thermal cycles are often taken as an indication for a first order phase transition. The effective 
potential allows to compute the values of p 2 also in the false vacuum and should hence repro- 
duce the hysteresis. In fig. |3] we show a thermal cycle at (5 = 8, A = 0.0001 on a 2 • 4 2 • 32 lattice. 
The data points, connected by a long dashed line to guide the eye, show a clear hysteresis. The 
solid line represents the values of p 2 in the absolut minimum obtained from the gauge invariant 
effective potential. The short dashed line indicate the values of p 2 in the second minimum. 
The figure demonstrates that the numerical data are very well described by the one-loop gauge 
invariant effective potential. The agreement gets worse for the high point at A = 0.0005. Here 
the transition points from the perturbatively evaluated potential and the simulations show some 
discrepancy, see table 0. 

In section [7] we will compute the latent heat. For this we need transition k's for L t = 2, 5. 
For the higher L t = 4, 5- values numerical simulations are very demanding as one would have 
to scale the other extensions of the lattice accordingly. Therefore we will resort there to the 
values of k c as obtained from the effective potential. For this purpose we performed a finite 
size scaling analysis of on various size lattices 

< = °y- v + «f , (21) 

where V = L x ■ L y ■ L z and L t is kept fixed. We computed for various V from the gauge 
invariant effective potential and fitted it to (PH). In all the fits performed, we found a value of 
v = 1.00(2). This again confirms the first order nature of the electroweak phase transition for 
Higgs masses below 50 GeV. The obtained results for for various L t , (3 and A are plotted 
in fig. T3 where we discuss the lines of constant physics. The numerical values in our four 
basic points are: 12 : «£° = 0.128290(1); 13 : «£° = 0.128082(1); h2 : «f = 0.128625(1); 
h3 : = 0.128273(1). 



4.2 Two-coupling method and transition points 

Provided hysteresis effects in thermal cycles are seen at a first order phase transition, the two- 
coupling method is useful for a precise determination of the position of the phase transition 
point. 

Let us consider the largest extension of the lattice to be the z-direction. In this direction the 
lattice is divided into two halves. The idea is to choose different coupling parameters in both 
halves, to enforce one part to stay in the symmetric phase and the other one in the Higgs phase. 
We assume that the z-direction is long enough for accomodating a pair of interfaces between 
the phases. (If the ^-direction is not too long, the appearance of more interface pairs has a 
negligible probability.) By warming up the configuration far enough in both directions from 
the transition point one can achieve that at the beginning of the simulation at the envisaged 
pair of parameters both phases and the interface pair are present. In this configuration the 
system can sensitively react to free energy differences by shifts of the interface positions in the 
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Table 2: Hopping parameter k c at the transition point, obtained with the two cou- 
pling method. As lattice size, the total size of both parts is given. For comparison, 
the estimates obtained on lattices with half the size from the one-loop invariant 
effective potential are also given. 



Lattice 


P 


A 


K c 


„1— loop 
K c 


2- 4 • 4 -128 


8.0 


0.0001 


0.12836(3) 


0.12840(5) 


2 -16 -16 -128 


8.0 


0.0001 


0.12830(5) 


0.12828(2) 


3 ■ 6 ■ 6 ■ 192 


8.15 


0.00011 


0.12813(2) 


0.12811(2) 


3 -12 -12 -192 


8.15 


0.00011 


0.12811(1) 


0.12811(1) 


2 -32 -32 -256 


8.0 


0.0005 


0.12887(1) 


0.12862(1) 


3 • 48 • 48 • 384 


8.15 


0.00051 


0.12852(2) 


0.12826(1) 



^-direction. If the configuration stays for many autocorrelation times in the mixed state the 
free energies of the two phases at the chosen parameters are such that this situation is stable 
against transitions to a unique phase. Thus the parameter sets of the two phases give a lower 
and an upper bound on the transition parameters. In fact, due to the additional free energy 
associated with the interfaces, the phase transition from a two-phase situation to a unique 
phase occurs somewhat earlier than the equality of the free energies of the two phases. This 
parameter shift goes, however, exponentially to zero for increasing lattice volumes. 

Of course, the two-coupling method for the determination of the transition point works only 
if one is able to tell one phase from the other. A possible way is to perform hysteresis runs 
on lattices with the same extensions as half of the lattice for the two-coupling method. The 
hysteresis plots are used for the distinction of the two phases at a given k. 

As stated in the introduction of this section, we are interested in the position of the phase 
transition in the hopping parameter at k = k c for fixed (3 and A. Thus the only parameter 
chosen differently in both parts of the lattice is the hopping parameter. At each pair of k the 
system was observed for at least 10 autocorrelation times. k c was defined as the mean value of 
the best lower and upper bounds. The best estimates for the transition point k c obtained with 
this method are given in table |2]. 



4.3 Multicanonical method and transition points 

The most precise way to determine the transition point has been provided by a combination 
of the multicanonical method and K-reweighting |J7||. In this method first an order parameter 
distribution is generated through a Monte Carlo simulation. This distribution is then extrap- 
olated to nearby re's to find the transition point. The action with the additional logarithmic 
term (Si og in eq. flTOD) and the link- variable (L v in eq. (|24|)) have been considered as order 
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lattice: 2 <g> 8 <g> 8 <8> 128 
= 8 k= 0.1283 A= 0.0001 




Figure 4: The distribution of the action density obtained from a constrained- 
multicanonical simulation on 2 ■ 8 2 ■ 128 lattice. Note that the left hand peak corre- 
sponds to the Higgs phase and the right hand one to the symmetric phase. 

parameters. The transition point is determined by the equal height signal: k c is the hopping 
parameter value for which the heights of the two peaks in the distribution are equal. The 
systematic uncertainty of k c can be estimated by comparing its values obtained by equal height 
signal and equal area signal for different order parameters. 

The multicanonical simulation has been performed at k = k v , close to the real transition 
point Qk v — k c \ = O(10~ 5 )). The approximate transition point k v was determined by the one- 
loop invariant effective potential. The distribution of an order parameter at a nearby k', with 
all other couplings fixed, can be obtained by attaching a weight 

W (L V ) = e mL ^'- Kv) (22) 

to the different configurations. 

The action density distribution at the low point for k v = 0.1283 on a lattice with Q = 
2-8 2 -128 sites is shown in fig. [|. The distribution in neighbouring points is obtained by eq. fl2"2"|). 
The two peaks are of equal height at k c = 0.128307 (fig. 0). As it can be seen, not only the 
heights are equal but also the widths of the two peaks are quite similar, thus the equal height 
condition for k c is roughly equivalent to the equal area condition. At the same time the flat 
regime between the peaks is almost constant. This means that in a multicanonical simulation 
the two phases can mix with each other with an arbitrary mixing ratio. The supression is due 
to the interfaces between the phases. Taking L v as an order parameter, the transition point 
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lattice: 2 <g> 8 <8> 8 <g> 128 
(3= 8 k= 0.128307 A= 0.0001 




S log 

Figure 5: The reweighted distribution of the action density obtained from the data 
of the previous figure at k — n c — 0.128307. 

corresponding to the equal height signal is given by a slightly higher k—0. 128308. This leads 
to a transition point at k c = 0.128307(2 + 1), where the first number in brackets denotes the 
systematic the second one the statistical error estimate. Similarly for Q = 2 • 4 2 • 128 and for 
D, = 2 • 4 2 • 64 the transition points are at k c = 0.128366(3 + 1) and k c = 0.128367(3 + 1), 
respectively. 

5 Masses and correlation lengths 

Important characteristic features of any statistical physical system are the correlation lengths. 
At zero temperature their inverses give the masses of the low lying particles. Particularly 
interesting are also the (inverse) correlation lengths on the two sides of the phase transition. 

5.1 Zero temperature masses 

The physical Higgs mass M H can be extracted from correlators of quantities as the site variable 

R x = ^Tr (<p+<p x ) = p 2 x , (23) 
or, using <p x = p x a x , the link variables 

= \ Tt i^t+fPx^x) , L V)Xlt = ^Tr (<f+ +(i U x ^(p x ) . (24) 
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lattice: 16 016 016 032 




5 10 15 

timeslice 



Figure 6: Correlation function in the Higgs boson channel in point /i2 [16/90]. The 
curve shown is the best fit for timeslices between 1 and 11 with Mu = 0.2663 and 
a constant factor /# = 4.081 10~ 2 . The %-square of this fit is x 2 = 0.54. The 
statistical error of the correlation function at distance 1 is about 0.5%, at distance 
11 about 6%. 



Table 3: The parameter values of numerical simulations for determining zero tem- 
perature masses and Wilson loops. 



index 


lattice 




A 


K 


sweeps 


12 


12 3 • 24 


8.00 


0.00010 


0.12830 


150000 


13 


18 3 • 36 


8.15 


0.00011 


0.12810 


125000 


h2[12] 


12 3 • 24 


8.00 


0.00050 


0.12890 


160000 


/i2[16/85] 


16 3 • 32 


8.00 


0.00050 


0.12885 


60000 


/i2[16/90] 


16 3 • 32 


8.00 


0.00050 


0.12890 


160000 


h3 


18 3 • 36 


8.15 


0.00051 


0.12852 


50000 
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Figure 7: Correlation function in the W-boson channel in point hi [16/90]. The 
curve shown is the best fit for timeslices between 3 and 15 with M w = 0.4522 and 
a constant factor f w = 1.946 10~ 5 . The x-square of this fit is x 2 — 1-30- The 
statistical error of the correlation function at distance 3 is about 0.2%, at distance 
15 about 26%. 



Table 4: The W-boson mass Mw and Higgs boson mass Mh in the points defined 
by the previous table. Their ratio is Rhw = Mh/Mw which leads to the Higgs 
boson mass in physical units given in the last column. 



index 


M w 


M H 


Rhw 


T c /M w 


M H (GeV) 


12 


1.059(24) 


0.236(7) 


0.222(12) 


0.472(11) 


18 


13 


0.718(3) 


0.144(4) 


0.201(5) 


0.464(2) 


16 


h2[12] 


0.427(8) 


0.262(9) 


0.614(32) 


1.171(22) 


49 


/i2[16/85] 


0.427(5) 


0.253(8) 


0.593(19) 


1.171(19) 


47 


/i2[16/90] 


0.453(3) 


0.266(5) 


0.587(12) 


1.104(7) 


47 


h3 


0.289(4) 


0.175(7) 


0.606(24) 


1.153(16) 


48 
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The W-boson mass M^r can be obtained similarly from the composite link fields (r, k — 1, 2, 3) 



Due to lattice symmetry only the diagonal correlators in (r, k) are non-zero and they can be 
averaged before extracting masses. 

Numerical simulations to determine the zero temperature masses were performed on 12 3 • 24 
and 16 3 • 32 lattices at the phase transition points of L t = 2 lattices. The former were scaled 
up to 18 3 ■ 36 at the corresponding L t = 3 points. Since the gauge-Higgs system is in the Higgs 
phase at low temperatures and in the symmetric phase at high temperatures, the T = points 
at the transition points are always in the Higgs phase. We did neither attempt to investigate 
the particle spectrum at T = in the symmetric phase nor did we determine the K-dependence 
of the masses in the Higgs phase in a wider range of k. These questions could be addressed in 
future studies. The only information on ^-dependence is obtained from two points on 16 3 • 32 
lattices at nearby K-values. The collection of parameters and lattices where zero temperature 
masses were determined is contained in table |3|. (In the same simulations the Wilson loops were 
also determined: see section ||.) 

Masses were extracted from the connected correlators by least square fits by a single cosh 
or cosh + constant. Constant contributions are possible in the Higgs channels but were most 
of the time negligible, because they are of the order exp(— M H L t ) and our lattices usually have 
large enough time extension L t . Typical examples in the /i2[16/90] point on 16 3 • 32 are shown 
by figures |6] and 0. 

Statistical errors on masses were always estimated by subdividing the data sample into 
subsamples. Performing the fits in subsamples gives estimates of standard deviations of fit 
parameters. This is particularly straightforward on the Quadrics Q16, because these lattices 
can be simulated on 8 nodes which is repeated 16 times with independent sequences of random 
numbers. The 16 parallel sets of statistically independent results can be used for the estimate 
of fit parameter errors. The error estimates in last digits are always given in parentheses. The 
obtained results are shown in table |j. 

The masses are in every case well determined. The correlators are always dominated by the 
lowest state. The statistical errors of the ratios of Higgs- and W-boson masses show that these 
two masses are not strongly correlated in the data samples. The values of T c /M w in table £| are 
obtained under the assumption that the points are at the finite temperature transition points 
of Lt = 2 and L t = 3 lattices, respectively. Therefore T c — 1/2 and T c — 1/3, respectively. 

Points /i2 [16/85] and /i2[16/90] are at slightly different k values: k = 0.12885 and k = 
0.12890, respectively, but otherwise identical. Their comparison shows that the ratio of Higgs- 
to W-boson masses remains almost constant in this K-range. The same holds also for the 
renormalized gauge coupling (see next section). The only significant difference is a small change 
in the overall scale shown by an about 5% decrease of the masses in lattice units. This feature 
is rather helpful because it makes fine tuning of the scalar hopping parameter k less relevant. 




(25) 
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Comparing points ft,2[12] with h2 [16/90] we see that the finite volume effects are not large. 
The only statistically significant deviation is seen in the W-boson mass, which is about 5% 
smaller on the 12 3 • 24 lattice than on 16 3 • 32 (for the change of the static potential see 
section ||). This implies that within our statistical errors the 16 3 • 32 lattice can probably be 
considered to represent the infinite volume limit. Nevertheless a systematic study of finite 
volume dependences in future studies is desirable. 

5.2 Correlation lengths at transition points 

An important effect of finite temperatures is the change in correlation lengths. This is particu- 
larly interesting at a first order phase transition, where correlation lengths stay finite and may 
eventually differ in the two metastable states. The inverse correlation lengths at finite temper- 
atures may also be called "temperature dependent masses" . In order to distinguish them from 
the real masses we shall denote them by small letters: mw for the inverse correlation length in 
the W-boson channel, tuh the one in the Higgs-boson channel. 

In our numerical simulations we make use of the fact that on a large enough lattice the 
metastability at the first order phase transition becomes so strong that in a finite amount of 
computer time the system stays in the phase it started from. At the first sight this seems to 
lead to mathematically ill-defined averages, because the transition speed of the lattice from one 
phase to the other does certainly depend on the particular updating algorithm. However, as it 
is well known (and is discussed in detail in section |j) the two phases can be well distinguished 
in the distribution of some order parameter which shows two well separated peaks on a large 
enough lattice. The averages in the two metastable states can be defined to belong to the two 
peaks in the order parameter distribution. The advantage of using large lattices is that one can 
sit right at the phase transition, hence no extrapolation is necessary (see also section ^|). This 
is more important than saving a little on lattice sizes and staying in an unclear situation with 
an underdeveloped two peak structure, when no clear separation of the two phases is possible. 

The lattice shapes we consider are typically of the type L t ■ L 2 xy ■ L z with L t = 2,3 the 
inverse temperature direction, L xy ^> L t the "small" spatial directions and L z 3> L xy the 
"long" spatial direction. The finite volume effects are mainly determined by L xy . The long 
direction L z is chosen to force the interfaces between the two phases to be built perpendicular 
to this direction. This is useful for the study of the interface tension (see section §). The 
parameters of our numerical simulations are collected in table |5|. The results are contained in 
table ||. 

The behaviour of the correlators in the Higgs phase (points with index a) is in general sim- 
ilar to the zero temperature points discussed in the previous subsection. There is, however, a 
substantial difference in the symmetric phase: in general the determination of the inverse corre- 
lation lengths is much more difficult. The main reason is that the correlators are not dominated 
by the contribution of the lowest state. There is obviously a dense spectrum contributing in 
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Figure 8: Correlation function in the Higgs boson channel in point I2[b] in the 
symmetric phase. The curve shown is the best fit for timeslices between 12 and 
28 with uifj = 0.2400 a constant factor /# = 4.287 10~ 3 and an additive constant 
ch = 2.974 10~ 5 . The %-square of this fit is \ 2 = 0.84. The statistical error of the 
correlation function at distance 12 is about 3%, at distance 28 about 20%. 



Table 5: The parameter values of numerical simulations for determining inverse 
correlation lengths at the phase transition. Indices [a] and [b] refer to simulations 
started "above", respectively, "below" the phase transition in /t-parameter, that is 
in the Higgs phase, respectively, symmetric phase. 



index 


lattice 


P 


A 


K 


a-sweeps 


6-sweeps 


I2[a,b] 


2 • 16 2 • 64 


8.00 


0.00010 


0.12830 


320000 


960000 


I3[a, b] 


3 • 24 2 • 96 


8.15 


0.00011 


0.12810 


144000 


80000 


h2[a,b] 


2 • 64 2 • 128 


8.00 


0.00050 


0.12887 


13000 


14000 


h3[a, b] 


3 • 96 2 • 192 


8.15 


0.00051 


0.12852 


5000 


5000 
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Table 6: The inverse correlation length in the W-channel mw and in the Higgs 
channel m H at the phase transition in the points defined by the previous table. 
Their ratio is raw = m H /raw- The physical scale is set by the last column, where 
the ratio of mw to the corresponding zero temperature mass Mw is given. 



index 


mw 


m H 




m w /T c 


m w /M w 


I2[a] 


0.880(15) 


0.130(3) 


0.148(6) 


1.76(3) 


0.83(3) 


I2[b] 


0.55(4) 


0.24(3) 


0.44(9) 


1.10(8) 


0.52(5) 


I3[a] 


0.569(8) 


0.0630(27) 


0.111(7) 


1.71(3) 


0.80(2) 


Z3[6] 


0.44(5) 


0.19(2) 


0.43(9) 


0.88(10) 


0.62(7) 


h2[a] 


0.218(6) 


0.0784(7) 


0.360(13) 


0.436(12) 


0.50(2) 


h2[b] 


0.34(6) 


0.101(12) 


0.30(9) 


0.68(12) 


0.80(16) 


h3[a] 


0.176(15) 


0.059(4) 


0.34(3) 


0.53(4) 


0.61(6) 


hZ[b\ 


0.34(7) 


0.0635(18) 


0.19(5) 


1.0(2) 


1.2(3) 
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Figure 9: Correlation function in the W-boson channel in point I2[b] in the sym- 
metric phase. The curve shown is the best fit for timeslices between 6 and 10 with 
m w = 0.5523 and a constant factor f w = 1.130 10" 5 . The ^-square of this fit is 
X 2 = 0.28. The statistical error of the correlation function at distance 6 is about 
4%, at distance 10 about 35%. 
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both channels. This feature is particularly pronounced in the W-channel. For this reason a 
rather high statistics was taken in the point with index I2[b] (see table |5]). Even in this case the 
lowest mass in the W-boson channel is still rather unprecise. For an illustration of the errors 
and fits see figures || and §. 

There is a clear signal for a discontinuity of m\y and m# between the two sides of the 
phase transition, whereby the most significant change is an increase of m# when going from 
the Higgs to the symmetric phase. The values of mw,H are rather small, therefore the lattice 
volumes may be too small, especially at the 12, 13 points. There the infinite volume values 
may still be somewhat different. From this point of view the situation is much better at the 
"high" points, where the volumes had to be chosen much larger due to the smaller value of the 
interface tension. 



6 Renormalized couplings 



In the SU(2) Higgs model there are two independent dimensionless renormalized couplings which 
are conventionally defined at zero temperature. Corresponding to the bare gauge coupling g and 
bare quartic coupling A one can introduce the renormalized couplings g R , respectively, Xr. One 
combination is conventionally fixed in the Higgs phase by the ratio of masses Rhw = Mh/Mw 
as 

R hw = —-r- ■ ( 26 ) 
9r 

This corresponds to the relations of the masses to the renormalized vacuum expectation value: 
M H = 8X R v R , respectively, = g R v R /4. For the correct normalization of Xr one has to 
remember that it corresponds to the bare coupling in perturbation theory Aq, which is related 
to A by 

^ (27) 

and Ao is often multiplied by 24 = 4!. 

Note that for fixed gauge coupling g R , because of the Weinberg-Linde bound on the Higgs 



boson mass f5B|| , Rhw an d X R defined by fl26|) have positive lower bounds. 

The mass ratio Rhw has been determined in the previous section (see table |J). Now we 
shall introduce the renormalized gauge coupling g R in a way which is convenient for numerical 
simulations. 



6.1 Renormalized gauge coupling 

The static potential at distance R can be obtained from the Wilson loops W(R,T) by 

V(R) = - Km hogW(R,T) . (28) 

1 — >oo ]_ 

We determined in the points defined by table |3] the on-axis Wilson loops of size R®T with 
1 < R < L s /2 and 1 < T < L t /2 on L? s ■ L t lattices. Every rectangular Wilson loop with two 
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Figure 10: The T-dependence of the Wilson loop ratio W(R,T — 1)/W(R, T) in 
the h2 [16/90] point at R = 4. The curve is the best fit with three exponentials. 



sides in time direction was included in the statistics. These Wilson loops were evaluated after 
transforming the gauge configuration to temporal gauge. The time dependence has been fitted 
by three exponentials in order to determine the large time asymptotics. We checked that by 
leaving out one exponential and going to larger time distances the lowest energy Vo does not 
change. In figure [10] we show a typical example for the behaviour of the logarithm of the ratio 
W(R,T — 1)/W(R,T). As it can been seen from the given fit parameters, the ground state 
and the excited states are separated by a large energy gap. The situation for all other points 
is quite similar. Due to this, there is no need for an optimization of the source to enhance the 
ground state signal. 

The static potential can be fitted well by a Yukawa-term with lattice corrections, as discussed 



in [pl |. It takes the form 

V(R) = -4 e " A/i? + C + D G(M, R, L s ) . (29) 
R 

The parameter M is the screening mass which is closely related but not exactly equal to 
the physical W-boson mass My/ determined from the correlation functions in section |5.1| . 
G(M, R, L B ) describes the difference between the continuum potential and the finite lattice 
version to lowest order. For on-axis R values it is given by 

G(M, R, L s ) = ^e~ MR - I(M, R, L.) , (30) 
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Figure 11: The continuum potential extracted (left picture) and the lattice cor- 
rection term (right picture) in points 12 and 13 as a function of distance. M is the 
screening mass. The dotted curve is the best fit for Z3. 



where 

47T e lkzR 
W L,)^£ M2 + 4ELsinW2) , (31) 

2imi 

hi = — — , Hi = 0, . . . , L s - 1 . (32) 

Lis 

In figure [ll| we show the continuum potential V cont (R) = V(R) — C — D G(M, R, L s ) and the 
correction term 5V(R) = D G(M, R, L s ) for the two points 12 and 13. One can see that on 
the L t = 3 lattice (point Z3) the lattice artifacts are at larger physical distances smaller but at 
R — 1,2 larger than on the L t = 2 lattice (point 12). 

The renormalized gauge coupling can be defined from the parameters of the static potential 
as g\ = 16ttA/3. In table [7] we collected the best fit parameters together with this global 
definition of g\. Another local definition of g 2 R) which is shown in the last column, will be 
described below. Obviously the renormalization effect on bare coupling g 2 = 0.5 is moderate. 
The nearby equality of A and D shows that the potential is dominated by the one vector boson 
exchange ||39|| . One can see also that finite volume effects are small (points ft,2[12] and h2 [16/90]) 
and that a small change in k (points ft,2 [16/85] and /i2[16/90]) leaves g\ almost constant. 

Another way to define the renormalized gauge coupling was discussed in We use a 

slightly different version which is similar to the procedure recently proposed for pure gauge 
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Table 7: Summary of the fit parameters for the static potential and the renormalized 
gauge coupling for both definitions. 



index 


A 


M 


D 


C 


9% = f ^ 


9 2 R (M- 1 ) 


12 


0.0337(3) 


1.113(16) 


0.0344(23) 


0.066453(6) 


0.564(6) 


0.563(6) 


13 


0.03340(6) 


0.730(4) 


0.0347(12) 


0.078005(3) 


0.5597(10) 


0.5611(23) 


h2[12] 


0.03434(7) 


0.4273(21) 


0.0352(8) 


0.090768(18) 


0.5754(11) 


0.5782(25) 


/i2[16/85] 


0.03440(6) 


0.4130(12) 


0.0365(8) 


0.091401(5) 


0.5763(10) 


0.5822(17) 


fc2[16/90] 


0.03431(3) 


0.4357(7) 


0.0373(4) 


0.090571(3) 


0.5749(5) 


0.5829(11) 


h3 


0.03373(8) 


0.2538(15) 


0.0366(9) 


0.095584(19) 


0.5651(13) 


0.570(7) 



theory [HJ. Given the potential, the renormalized coupling at distance Ri is defined by 

n 2 (R) = 1 ^ V(R)-V(R-l) m] 
9rW- 3 I{M ,R-1,L S )-I(M,R,L S ) ' {66) 

This can be interpolated to a physical distance, for instance r = M _1 . Ri is defined by the 
force as the solution of the equation 



}_ e -MRi 



— + M 
Ri 



I(M,R-1,L S )-I(M,R,L S ) . (34) 

ttj L it j 

The value of M in eq. (|33| ) can be taken from the fit ( |29| ) restricted to large distances. The reason 
of introducing a local definition of g R is that the short distance potential is not purely Yukawa 
like. Specifying the distance is therefore in principle important. However, our resolution is up 
to now not fine enough to see the logarithmic dependence on distance. 

The time dependence of the Wilson loops as well as the distance dependence of the static 
potential were described very well by the ansatze, i.e. y 2 was of order one in both cases. The 
statistical errors were determined from independent subsamples (see previous section). The 
error for the renormalized gauge coupling squared g^(M _1 ) quoted in table [7| is the sum of 
the statistical errors for potential and mass, and the systematical error. The interpolation to 
the distance M _1 was done by a linear approximation of the two neighbouring points. The 
systematical errors are estimated from the change arising by the inclusion of a third point. 
The given errors are usually dominated by the statistical errors of M. Note that the value of 
^(M -1 ) in point 12 is corrected as compared to [2C]. 



6.2 Renormalization group trajectories 

The continuum limit of quantum field theories on the lattice is taken along renormalization 
group trajectories, also called lines of constant physics (LCP's). Going to the continuum limit 
along such lines in bare parameter space the same physical theory is reproduced with increasing 
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0.2 0.4 0.6 0.8 1 
•r-iog(L t /2) 

Figure 12: The positions of the phase transition in k along two lines of constant 
physics, which start at the phase transition points of an L t = 2 lattice at (3 — 
8, A = 0.0001 (lower curve) and (3 — 8, X = 0.0005 (upper curve). The points 
with 2 < L t < 5 are determined from the one-loop invariant effective potential. The 
curves are polynomial interpolations. 



precision. This means that the physical results of numerical simulations should be the same, 
apart from "scale breaking lattice artifacts". In the SU(2) Higgs model, in order to define 
the LCP's, one has to keep fixed the values of two independent renormalized couplings, say 
the above discussed renormalized quartic (Xr) and gauge (g^) couplings. The third parameter 
characterizing the points along LCP's can be chosen as 



T EE logM^ 1 



(35) 



(Remember that in this paper the lattice spacing is set to a = 1, therefore Mjy is measured 
here in lattice units.) 

Since our numerical simulations are performed in the weak coupling region of parameter 
space, the change of the bare couplings g 2 = 4/(3 and A defined by eq. ( p7|) along LCP's 
can be well approximated by the solutions of the one-loop perturbative renormalization group 
equations: 



dg\r) 
dr 



16tt 2 
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g* + O(XlXlg 2 ,X g\g 



dXo(r) 
dr 



16tt 2 



9 



96A ° + 32^ ~ 9X ° 9 + ° (A °' 09 ' X ° 9 ' 9 



(36) 
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The integration of these equations can be started at the transition points of the L t = 2 lattices. 
At a distance At = log(L t /2) (L t = 3,4,5, . . .) in parameter r one obtains the corresponding 
values of (g 2 , Xq) for the transition points of lattices with temporal extension L t . 

In order to have a perturbative prediction also for the change of the third bare parameter k, 
one can make use of the one-loop perturbative invariant effective potential. The infinite volume 
predictions for our "low" and "high" values of the quartic coupling in the range 2 < L t < 5 
are shown by fig. [12|. As we discussed in section f|, the measured values of k c are very well 
reproduced by the one-loop invariant effective potential for L t = 2, 3 at our "low" value of the 
quartic coupling. At the "high" quartic coupling there is a parallel shift between prediction 
and measurements: the prediction is too low by Ak c = 0.00025(1). As a consequence, the 
derivatives of k along these LCP's can be well determined: 

dn 12 dn h2 
—— = -0.000685(40 + 69) , — — = -0.00120(9 + 5) , 
Or or 



cW 3 . dn 



h3 



-0.000367(17 + 4) , — — = -0.000588(35 + 6) . (37) 
or or 

The errors given here include an estimate of the systematic error by comparing different poly- 
nomial interpolations of the points 2 < L t < 5 (first entry in parentheses) plus the error coming 
from the uncertainties of k c obtained in the numerical simulations (second entry). 

Although the estimates of the LCP's by one-loop perturbation theory are very useful for 
guiding the numerical simulations, finally one has to confront them with numerical simulation 
data. Having this in mind we have chosen the L t = 3 points with index 13 and h3 on the 
trajectories of the solutions of eq. ( |36"D which start at 12 and h2, respectively. (For j3 and 
A we kept only the first few digits which are sufficient within our typical statistical errors.) 
Comparing the values of Rhw an d g 2 R in the points 12 and /3, respectively, h2 and h3 (see 
tables f| and 0) one can conclude that the one-loop predictions from eq. ([36|) are very good: the 
measured values of both renormalized couplings in corresponding points are very close to each 
other, in most cases equal within errors. 



7 Latent heat 

An important characteristic feature of first order phase transitions is the latent heat, i.e. the 
discontinuity of the energy density e. The pressure P is continuous, therefore Ae can be 
obtained from the discontinuity of 5 = e/3 — P. This latter quantity is related to the trace of 
the energy-momentum tensor and is theoretically simple, because it can be obtained from the 
derivative of the action density with respect to the lattice spacing a. This, in turn, implies that 
it is closely related to the LCP's for the continuum limit discussed in the previous section. 



As it has been shown in |20|, the latent heat Ae in the SU(2) Higgs model is given by 



% = 14 (| • - | • AO. - f • 6AP pl ) . (38) 
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Table 8: The obtained values of global averages "above" and "below" phase tran- 
sition points, i.e. in the Higgs, respectively, symmetric phase. 



index 


^pl 




T 


r 


^XX 


c 

*->x 


I2[a\ 


0.085487(5) 


25.184(12) 


0.91666(4) 


22.720(12) 


652(2) 


6.0333(3) 


I2[b] 


0.096053(1) 


2.88354(17) 


0.27614(3) 


0.86267(16) 


7.5(1) 


6.60939(2) 


I3[a] 


0.09009(12) 


11.548(24) 


0.8284(3) 


9.349(24) 


138.4(5) 


6.3878(6) 


I3[b] 


0.094438(1) 


2.6585(3) 


0.21724(7) 


0.6442(3) 


6.2802(17) 


6.61696(3) 


h2[a] 


0.095391(6) 


3.988(11) 


0.4629(15) 


1.947(11) 


16.17(9) 


6.5677(3) 


h2[b] 


0.096000(1) 


2.9761(15) 


0.2964(3) 


0.9568(15) 


8.302(10) 


6.60189(6) 


h3 [a] 


0.094183(2) 


3.133(3) 


0.3243(6) 


1.113(3) 


9.31(2) 


6.59861(9) 


h3[b] 


0.094421(1) 


2.7006(11) 


0.2290(3) 


0.6896(11) 


6.517(7) 


6.61214(4) 



Here, besides the notations in (H), we used 

Qx = (pI - 1) 2 , 



pi 



(39) 



The sign of the discontinuities like Ae etc. will be defined in such a way that they are differences 
of values in the symmetric phase minus Higgs phase. 

The obtained average values of the global quantities in eqs. (|23|) , ([24]), (|39|) and of the 
lattice action per point 

S x = 6f3P pl + R X + XQ X - %kL V w (40) 

are shown in table |8|. The indices of points where the numerical simulations were performed 
are defined in table |^. 

As it has been discussed above (see e.g. section |5]^), the small errors in table [8] could 
be achieved by avoiding errors of extrapolation. This is made possible by performing the 
simulations on large enough lattices, where the strong metastability of phases can be exploited. 
In this way the contribution of the statistical errors in table |8| to the errors of Ae is negligibly 
small compared to the errors of the derivatives. In particular, in our case the errors in eq. (|4T|) 
are dominated by the errors coming from eq. (|37|). 

Besides the discontinuities of global quantities, for Ae in eq. ( |3"8"D the derivatives of bare 
parameters along the LCP's are needed. These have been discussed in the previous section. 
dp/dr and dX/dr can be obtained to a good approximation from eq. (|36|) . The values of dn/dr 
have also been evaluated from the one-loop invariant effective potential and from numerical 
simulation data. They are given by eq. (37). Inserting them into eq. (38) together with the 
numbers from table 



one obtains 

12 




1.81(29) 




0.132(17) 
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(^j =1.57(13), ^ =0.122(9). (41) 

We see from here that within errors the results on L t = 2 and L t = 3 lattices coincide. 

Comparing the contributions of different terms to Ae in fl38|) one sees that the term with 
L v dominates. This is in accordance with the observation that the lattice spacing is mainly 
determined by the hopping parameter k (see sections f| and ||). 



8 Interface tension 

At the transition point of the electroweak phase transition mixed states can appear, where 
different bulk phases are separated by interfaces. The interface tension, a, is the free energy 
density associated with these interfaces. The dynamics of a first order phase transition is to 
a large extent determined by the interface tension. In particular, the nucleation of bubbles or 



droplets is essentially influenced by it [[41], [42j1 - In perturbation theory the determination of a 
is a task which is far from being trivial. The main problem is the bad infrared behaviour of 
the finite temperature effective action in the S'?7(2)-Higgs model. Therefore, non-perturbative 
lattice simulations are very useful. 

In this section we use two different methods, namely, the two-coupling method and Binder's 
histogram method combined with multicanonical updating, to determine a. In both cases we 
use elongated lattices. This geometric choice of the lattices results in interfaces which are 
perpendicular to the long direction. Other configurations (e.g. multiple walls, large spherical 
bubbles) are highly suppressed. The inverse temperature will be restricted in this section to 
L t = 2. 



8.1 Two-coupling method and interface tension 

One way to determine the interface tension in the SU(2) Higgs model is based on a modified 
version of the Potvin-Rebbi two-coupling method Similarly to subsection [4.2| , a lattice is 
considered with a long extension in the ^-direction. The long direction with periodic boundary 
conditions is divided into two equal halves with different couplings, in such a way that an 
interface pair is created near the hypersurfaces where the couplings are changing. 

Since we have three bare parameters, in principle one can choose any of them (or some 
combination) to be different in the two halves of the lattice. A physically good choice would 
be to change the temperature only. One can move, for instance, along an LCP (see section |6.2| ) 
when only the lattice spacing a is changing. Assuming that the volume is very large, the only 
relevant change is in the temperature T = l/(aL t ). A simpler way, which comes close to this, 
is to change only k. The reason is that in the latent heat the contribution of the y)-link L v 
dominates which is conjugate to the hopping parameter k. 

Therefore we decided to split the value of k. On an L t -L x -L y -L z lattice with L z ^> L x ^ y j, and 
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in our case L t = 2, L x = L y = L xy , the two halves in the z-direction have hopping parameters 
K± t 2 which are slightly below and above the transition point: K\ < k c < k 2 . 

The interface tension is defined by the "free energy" (more precisely "free action") 

F = S - ttd(s) = {lf(s) , (42) 

where S is the lattice action ([I]), Q = L x L y L z L t is the number of lattice points and d(s) is the 
spectral density of states as a function of the action density s = S/ft: 



e nd( s ) 



= JldU][d V ]s(s-^) . (43) 



This implies that the probability distribution of s is w(s) oc exp[— Qf(s)}. The interface tension 
between the states with K\ and k 2 is given in lattice units by 



o" = (2L x L y L t ) 1 



1 1 

k 2 ) - 2 F ( Kl > Kl ) ~ 2 F ( K2 > 



' {[F(ki, k 2 ) - F(ki, ki)] — [F(k 2 , k 2 ) — F(k u k 2 )]} . (44) 



AL x L y L t 



Let us write the action on the lattice with two halves as 

S = S - /tiS 1 ! - K 2 S 2 , (45) 
where So is the piece independent from k. Then in the thermodynamical limit we have for 

^fe^ —<*>—. («) 

if (. . .) K1)K2 denotes expectation values for a lattice with two halves at Ki )2 . Therefore we obtain 
from ( p|) 

a(/ci,K 2 ) = {^L x LyL t y l \ j dK{S x ) K>K2 - / ck(S 2 ) Kl)K \ . (47) 

Here the notation cr(/«i,/« 2 ) emphasizes the dependence on « lj2 . Finally, since S i)2 is given by 
the yj-link L v , we have 

<t(ki, k 2 ) = L z \ f 2 dKL$\K, k 2 ) - I 2 dKL$\K U re)} , (48) 

where iA 1 ' 2 -^/?, re') denote the expectation values of averaged in the two halves, if the hopping 
parameters are k and re', respectively. 

In order to obtain the physically interesting interface tension one has, of course, to perform 
the non-commutative limits L Xl y, z — > oo and Are = re 2 — re c = re c — rei — > 0. For a given 
lattice extension Ak cannot be arbitrarily small, because if the difference in free energy density 
becomes small tunneling into the other phase can occur and the interfaces disappear. The 
presence of the interfaces can, however, be monitored to ensure the applicability of (§3). 
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For small Ak the integrals in (|4q) can be well approximated by the average of the integrand 
at the two end points. Besides, for equal arguments we obviously have L^'(k, k) = Ly(n, k). 
This gives 



cr ~ L 7 Ak 



L^\k 1 ,k 2 )-L^\k 1 ,k 1 ) + L^(k 2i k 2 )-L^(k 1 ,k 2 )\ . (49) 



We are interested in the limit Ak — > 0. Therefore in the square bracket only contributions of 
order 1/Ak matter. Such contributions cannot come from terms with equal hopping parameters 
in the two halves. Therefore let us introduce the parametrizations 



L^\k u k 2 ) = - h h + ai(«! - k c ) + 0{k x - k c ) 2 

K\ K c 



L^(k u « 2 ) = — + b 2 + a 2 {n 2 - k c ) + 0(k 2 - k c ) 2 . (50) 

K 2 - K c 

Then for Ak — > a finite volume estimator of the interface tension is 

a = L z ( Cl + c 2 ). (51) 

In view of this formula the question arises about the origin of the contributions of order 
0(1/ Ak). For small Ak the difference in free energies of the two phases is of the order O(Ak). 
Therefore the interfaces are not fixed at the hypersurfaces where k changes, but are penetrating 
into the neighbouring regions with constant k. If the k change is at z = zq, the probability of 
the interface position at z > z is oc exp[— const. Ak(z — Zq)\. Integrating over z gives 0(1/ Ak). 
In fact, as the numerical simulation data show, the distributions of ^-slices of as a function 
of z can be well fitted by exponentials. 

The above formulas also show that the extrapolation to Ak = is delicate. The intervals 
in k 12 for the fits in eq. (|5T)D have to be carefully chosen. These forms are, in fact, replacing 
(fHf) in the region where <j(ki,k 2 ) depends linearly on K\ )2 . For given lattice extensions one 
cannot take intervals too close to k c . This is a kind of round-off effects, which usually appear at 
first order phase transitions in finite volumes but get smaller for increasing lattice sizes. First, 



according to eq. ([51]) one has to choose L z large enough in order that the two interfaces do 
not interact with each other. In this case the left hand side becomes independent from L z . 
Second, for increasing L x , L y the free energy differences grow proportional to (L x L v Ak). The 
probability of penetration into the regions with constant k is exponentially decreasing. This 
allows to consider small Ak's of the order of 0(l/(L x L y )). 

In the numerical simulations with two k's the lattice sizes extended up to 2 • 16 2 • 128 in 
the 12 point and 2 ■ 32 2 ■ 256 in the h2 point. The numerical simulation data on these lattices 
are collected in table [| respectively, table |10[ Fits with the parameters a, b, c in eq. (|50|) 
were performed, omitting quadratic and higher order terms in Ak. For the transition points 
k c = 0.12830 and k c = 0.12887 were assumed in points 12 and h2, respectively. The data entries 
in the tables with italic were not taken into account in the fits, in order to obtain acceptable 
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Table 9: The average values of L^' 2 ^ obtained in two-n simulations with k = Kip 
in point 12 on 2 • 16 2 • 128 lattice. 



Ki 


«2 




4 2) 


0.12790 


0.12870 


2.6832(44) 


53.422(12) 


0.12795 


0.12865 


2.7148(54) 


50.088(17) 


0.12800 


0.12860 


2.7372(65) 


46.630(16) 


0.12805 


0.12855 


2.7731(84) 


43.052(15) 


0.12810 


0.12850 


2.8323(87) 


39.300(17) 



Table 10: The average values of L^' 2 ^ obtained in two-K simulations with k = k 1j2 
in point h2 on 2 • 32 2 • 256 lattice. 



Ki 


«2 




4 2) 


0.12880 


0.12894 


0.9948(20) 


2.9325(43) 


0.12881 


0.12893 


0.9991(29) 


2.7939(65) 


0.12882 


0.12892 


1.0173(24) 


2.6688(54) 


0.12883 


0.12891 


1.0269(54) 


2.5217(86) 


0.12884 


0.12890 


1.0531(85) 


2.377(10) 


0.12885 


0.12889 


1.0873(92) 


2.237(13) 


0.12886 


0.12888 


1.043(11) 


1.998(28) 
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X 2 's. The other points can be well fitted and the x 2 's turn out to be of the order of the number 
of degrees of freedom. Inserting the results into eq. ( CT ) and assuming T c = 1/2 we obtain 



l^j =0.84(16), l^j = 0.008(2). (52) 

The errors here also contain some subjective estimates of the systematic errors coming from 
the choice of fit intervals. The statistical part of the errors was estimated by repeating the fits 
with normally distributed input data. 

Concerning the volume dependence of a we also collected data in point 12 on 2 ■ 8 2 • 128 
lattice. Within statistical errors of about 10% no deviation from the result given in (|52"D could 
be detected. 

8.2 Multicanonical method and interface tension 

In this subsection we present our results for the interface tension obtained by the histogram 
method |44[] at the low point, L t = 2. 

At k c the probability distribution of an order parameter (e.g. action density si og = Si og /Q 
defined by (|T0"D, or link variable L v in (p4|)) develops two peaks. They correspond to pure phases 
and the suppressed configurations between the peaks are dominated by mixed states where the 
phases are separated by interfaces. Defining n c by the equal height signal the suppression at 
infinite volume is given by the interface tension 

o-oo = lim cr n , o-Q, = nT \ T lo S , ( 53 ) 

\l~ >oo LL x L y L t Pmin 

where p max corresponds to the heights of the peaks and p m in to the minimum in between. 
Recently intensive studies have been carried out in order to understand the finite size cor- 



12 / „ \ h2 



rections to a [45, 46, 47, 45, 45]. We follow [48] and for our simulations at the "low" value of the 
quartic coupling A = 0.0001 use elongated lattices with extensions Q = 2 • 4 2 • 64, Q = 2 • 4 2 • 128 
and Q = 2 • 8 2 • 128. In this case the finite size corrections are particularly simple 

13 1 
o-oo = cn + 75-jr (c+ -rlogL,- -logL xy ) , (54) 

L xy L t 4 I 

where L x = L y = L xy , L z is the longest extension and c is an unknown constant. The obvi- 
ous advantage of this choice is that practically all mixed configurations contain two surfaces 
perpendicular to the z direction. More than two surfaces or surfaces in any other directions 
are suppressed by many orders of magnitude. (Corrections beyond the gaussian approximation 
49| are negligible in our case.) 

The basic picture behind the above formulas assumes that the two interfaces are thin and 



do not interact. Fig. [13| is a typical picture of the average L v values for different z-slices in a 
mixed configuration. As it can be seen, in the case of Q = 2 • 4 2 • 64 the interfaces are very 
close to each other and they cannot be considered as thin (the spatial width in lattice units 
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lattice: 2 <g> 4 <g> 4 <g> 64 
= 8.0 /c=0. 12837 A = 0.0001 s log =2. 80638 
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Figure 13: The average values of on z-slices in a configuration of 2 • 4 2 • 64 
lattice. The dotted lines indicate the positions of maxima of the L 9 distributions in 
the bulk phases. 
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25 
20 
15 

a. 

10 

5 





50 100 

z 

Figure 14: The average values of on ^-slices in a configuration of 2 • 8 2 • 128 
lattice. 
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is approximately 30). In this case there is no flat region between the two peaks of the Si og 
distribution. The situation is much better for the even longer lattices (e.g. tt = 2 • 8 2 • 128, 
see fig. |4]). In the latter case the distribution of the action density si og at shown in 

fig. [|. It is remarkably flat between the peaks. This is another sign that the interfaces are well 
separated, thus mixed configurations can appear with any bulk phase ratio. 

In order to suppress statistical noise we have used a third order polynomial fit to the 
histograms near the extrema and calculated the extremum values from the fits. (However, due 
to our good statistics the simple use of the histogram-bins has resulted in almost the same 
results.) Since in case of multicanonical simulations one lattice has been attached to each node, 
we had 128 independent simulations. Our error estimates have been obtained by inspecting 
these independent samples. 

The finite volume interface tensions turn out to be ctq/T c 3 = 0.63(3) for Q = 2 • 4 2 • 128 and 
o^ijTl = 0.80(2) for Q = 2 • 8 2 • 128, respectively. Combining these two values one gets from 

~ 0.83(4) . (55) 

This result shows that the infinite volume limit is not far away from our largest volume result 
and it is in very good agreement with the results of the previously discussed two-ft method in 




9 Summary and discussion 

The main conclusion of this paper and ref. [^0J is that the electroweak phase transition in 
the SU(2) Higgs model is of first order for Higgs boson masses below 50 GeV . The strongest 
indication for this is the two-peak structure of order parameter distributions (see e.g. fig. |5|) 
showing a strongly suppressed flat region between the peaks which corresponds to mixed phases. 
For Mu — 50 GeV the phase transition became rather weakly first order. Still we could simulate 
large enough lattices to reach long living metastability. In fact, we used this phenomenon to 
study the properties of both metastable phases separately (see e.g. section |5\2| ). Investigations 
of the correlation lengths in both phases showed that they behave discontinuously at the phase 
transition. The strength of the first order phase transition can be characterized by the latent 
heat (section 0) and by the interface tension (section ||). Both these quantities play a central role 
in the dynamics of first order phase transitions and can be determined in numerical simulations. 
Our results show that at small Higgs boson masses (Mg ~ 20 GeV) the transition is very strong, 
with a large latent heat (eq. (f4T|)) and interface tension (eqs. (|52|), (|55l)). At Higgs boson masses 
near M# ~ 50 GeV, however, the phase transition becomes weak: the dimensionless latent heat 
decreases by about one order of magnitude, the dimensionless interface tension by about two 
orders of magnitude (see the same equations). The rate of decrease is qualitatively the same as 
given by two- loop resummed perturbation theory ||. For this we refer to the figures of ref. [^(J 
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Figure 15: The ratio of the critical temperature and the T = Higgs mass M# 
as a function of M H . The two curves are the results from the two-loop continuum 
calculation for our renormalized couplings and Mw = 80 GeV. The numerical 
results are the squares with error-bars for M H ~ 18 and 49 GeV. The smaller value 
in both cases corresponds to L t = 2 and the larger one to L t — 3. The triangles 
represent the results obtained by the lattice version of the gauge invariant potential 
which also give growing T C /M H values. This growth is shown by the insert in the 
upper right corner. 
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which are very similar after the inclusion of the new results here. 

Since this is our first large scale numerical simulation of the electroweak phase transition, 
we could not yet achieve a full quantitative control of the errors. This definitely requires further 
essential work. Nevertheless, our lattices are large and the statistical errors in most cases small. 
In the present paper first attempts were also made to estimate the systematic errors due to 
finite lattice spacing and finite volume. The comparison of the results on L t = 2 and L t = 3 
lattices shows remarkable consistency at the 10 % level. This means that lattice artifacts are 
not large. In particular, the latent heat shows a remarkably good scaling between L t = 2 and 
L t = 3 lattices. The apparent general tendency of small scaling violations already at such small 
temporal lattice extensions is consistent with the expected dominance of the lowest Matsubara 
modes motivating the dimensional reduction [L6|, [IT]] . 

Comparing the results for T c from the lattice simulations and those from perturbation theory, 
the values obtained from the simulations are always smaller. However, the proper comparison 
should be done not with continuum but with lattice perturbation theory, for relatively small 
LtS. We have carried out this analysis (see fig. [15]). In order to determine the critical points 
and masses the gauge invariant effective potential of section |3| has been used (again, the role 
of the symmetric phase and effects due to higher order corrections were neglected). One can 
see that for m# = 18 GeV the perturbative result is in complete agreement with the lattice 
one. A scaling violation can be observed and it goes approximately like 1/L 2 oc a 2 . (See the 
inlet of fig. [15|, where L t « 5 — 7 practically reaches 2.47, indicated as a dashed line. This 
value corresponds to L t — > oo, which can be also obtained by the use of the one-loop effective 
potential without high-temperature expansion.) In the case of ran ~ 49 GeV the situation 
is somewhat different. The perturbative values are larger than those obtained by the lattice 
simulations. The scaling violation is much smaller and within the errors it is reproduced by 
the perturbative treatment. 

Up to now we did only a few checks of finite volume effects. These showed that the results 
are only changing at the few percent level. Further systematic tests have to be done in the 
future. Nevertheless, the lattices are large enough to support strong metastability of the two 
phases. Therefore the conclusions concerning the strength of the first order phase transitions 
are firm. 

In general, numerical simulations of the electroweak phase transition in the four dimensional 
SU(2) Higgs model in the range of Higgs boson masses below Mh — 50 GeV turned out to 
be feasible and powerful. For these relatively light Higgs masses there is no need to go to the 
three dimensional reduced model. The situation could, however, be different for heavier Higgs 
bosons. In fact, recent numerical studies of the reduced SU(2) Higgs model are concentrating 



on M H ~ M w 50 



In future numerical simulations in the four dimensional model one has to consider also 
L t = 4 lattices for confirming the smallness of scaling violations. A systematic study of finite 
volume effects is desirable. Of course, an extension towards heavier Higgs bosons would be very 
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interesting. 
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A Appendix: multicanonical heatbath algorithms 



Here we describe the modifications for the heatbath algorithms in multicanonical simulations. 
For sake of simplicity we first explain the general strategy in case of the gauge field and then 
the realization of the algorithms for the fields U and (p. We use the notations introduced in 
section |2l|. In general, the indices old and new will denote the value before the update and the 
proposal, respectively. 

We distinguish two cases. If all possible values of the new action remain in the interval 
h = (S k , S k+1 ], one has to replace S by (1 + (3k) S and the generation of the distribution is as 
usual (see section [2.1|) . If S^™ crosses an interval boundary, the two distributions which are 



different for the two intervals have to be matched. This is done in such a way that for the interval 
Ik the distribution becomes Wk ~ exp[—(l + (3k)S] and for Ik+i it is Wk+i ~ exp[— (1 + (3k+i)S}. 
In addition the distribution for the combined interval Ik and Ik+i have to be continuous. The 
idea is to construct a distribution which is a hull of Wk and Wk+i- To obtain the desired 
distribution in each interval one has to correct the hull distribution by imposing an accept- 
reject step. The correction factor (denoted by 0) has to be positive and not larger than one. 
Clearly if the hull distribution is close to Wk and Wk+i the correction factor is about one and 
the acceptance rate is high. 

A.l Gauge field 

Let's first discuss the gauge updating. For the standard heatbath algorithm the main task is 
to generate a distribution 

P(a )da ~ (1 — a 2 )2 exp(l;a )da (—1 < a < 1, £ > 0) . (56) 

This is done by the algorithm described in ||26|| . In the multicanonical case the first step is to 
check, as discussed above, whether the global action SJ^° can cross an interval boundary. This 
depends on a through S?™{a% ew ) = Sf Q d + £« d - a% ew ). If the minimal action 5^(1) and 
the maximal action SJ^l—l) are in the same interval Ik one replaces £ by £(1 + (3k) in fl5T)| ) and 
proceeds in the usual way. In case that the action crosses an interval boundary S h+1 , which 
defines the corresponding a b Q Und through S k+1 = (a b Q Und ) , one constructs the distribution 
in the following way: 

(a) (3k < Pk+l 

Generate the hull distribution 

P(a )da ~ (1 — a 2 )2 exp[£(l + (3k)ao]da . (57) 

If a > a b ound (Sjffi e h) no correction factor is needed. For a < a b ound (S?™ e h+i) 
one has to impose 

6 = exp[£(/3 fc+1 - (3 k )(a - a bound )] (58) 
as an additional accept-reject factor (see figure jlB]) . 
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Figure 16: The figure illustrates the different cases for gauge field updating. In 
case Pk < Pk+i the desired multicanonical distribution is represented by the dashed 
line, in case Pk > Pk+i by the continuous line. In the first case one generates the 
distribution for the whole interval — 1 < Oq < 1, the continuation of which is 
the solid line for ao < a b ) ound , and corrects only if ao < a b Q Und . In the second case 
one generates a hull distribution (dotted line) proportional to in the whole 

interval. 



(b) p k > (3 k+l 

Generate the hull distribution 



P(a )da ~ (1 - a ) 2 exp[£(l + P k+ i)a ]da 



and correct for the following factor (see figure jig): 

exp^Pk - Pk+x)(ao - 1)] a > a bound (S?™ G 4) , 

„ ^- „bound I cnew r- j "> 

a < a \b log fc i k +i, 



e 



exp[aPk-Pk + i)(a b ound -l)] 



(59) 



(60) 



To avoid that the action Sf^ crosses more than one interval boundary, the intervals must be 
large enough. 

In practice the values of P^s are small (O(10~ 2 — 10~ 3 )). Therefore the two curves are closer 
together than in figure [16| and the acceptance rate is high. Also the values of £ are higher, i.e. 
there is a sharp peak near ao = 1. 
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A. 2 Scalar field 



The algorithm which is used for ^-updating was described in section |2J]. We write the action 
of eq. (Q) as a sum of the quadratic and quartic part and a constant 

S(ip x ) = S( 2 )((Px) + S(4)(ip x ) + const.' . (61) 

For the multicanonical simulations one has to generate the distribution 

P{(p x )d 4 ip x ~ exp[-S(ip x ) - P k Si og - a k ]d 4 ip x (62) 

for the four components of (p x . The first step is to generate the distribution 

P(<p x )d 4 <p x ~ exp[-(l + f3min)S(2)(¥x)}d 4 (p x (63) 

as a hull analytically by replacing S(2)(<p x ) — > (1 + ^min)^(2)(fx)- Pmin denotes the minimal fa. 
Let Ik denote the interval in which the minimum Sf of the action can be found. If Sf™ is in 
this interval Ik, the correction factor is 

0i = exp[(p min - Pk)S {2 )(Px) + 3 fa logGO - (1 + fa)S(4)(<Px) - c] (64) 

where the meaning of the constant c is described below. In case of crossing an interval border, 
i.e. SJZjj" = S — 3 J2 X l°g(Px) > S k+1 , one needs an additional factor 

9 2 = exp{[Sr o e g w - S k+1 ](fa - fa +1 )} . (65) 

The total correction factor is in this case 0i<92- This procedure generates the action distribution 
in Ik and in Ik+\ correctly. Due to the definition of J^, cannot be in Ij (J < k). Since our 
intervals are large enough the probability that SJ 1 ^" is in Ij (J > k + 1) is negligible. In practice 
it is quite difficult to find Sf og , but a slightly smaller Si og < Sf og can be used, which could lead 
to a somewhat smaller acceptance rate. In addition it is necessary that the whole correction 
factor (which includes all terms except the analytically generated quadratic term) has to be 
smaller or equal one. To be sure that no problem arises we subtracted a small positive extra 
term c in the exponent of Q\. This term reduces the acceptance. It must be tuned in such a 
way that the correction factor is never greater than one and the acceptance rate is high. 

Of course, there is always only one accept-reject step in both updating algorithms. This 
means that the analytically generated distribution will be corrected with the product of all 
factors. 
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